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ABSTRACT 

Context. After the termination shock (TS) crossing, the Voyager 2 spacecraft has been observing strong variations of the magnetic 
field and solar wind parameters in the heliosheath. Anomalous cosmic rays, electrons, and galactic cosmic rays present strong 
intensity fluctuations. Several works suggested that the fluctuations might be attributed to spatial variations within the heliosheath. 
Additionally, the variability of the solar wind in this region is caused by different temporal events that occur near the Sun and propagate 
to the outer heliosphere. 

Aims. To understand the spatial and temporal effects in the heliosheath, it is important to study these effects separately. In this 
work we explore the role of shocks as one type of temporal effects in the dynamics of the heliosheath. Although currently plasma in 
the heliosheath is dominated by solar minima conditions, with increasing solar cycle shocks associated with transients will play an 
important role. 

Methods. We used a 3D MHD multi-fluid model of the interaction between the solar wind and the local interstellar medium to study 
the propagation of a pair of forward-reverse shocks in the supersonic solar wind, interaction with the TS, and propagation to the 
heliosheath. 

Results. We found that in the supersonic solar wind the interaction region between the shocks expands, the shocks weaken and 
decelerate. The fluctuation amplitudes of the plasma parameters vary with heliocentric distance. The interaction of the pair of shocks 
with the TS creates a variety of new waves and discontinuities in the heliosheath, which produce a highly variable solar wind flow. 
The collision of the forward shock with the heliopause causes a reflection of fast magnetosonic waves inside the heliosheath. 

Key words, magnetohydrodynamics (MHD) - shock waves - waves: heliosphere - solar wind - Sun 
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The heliosheath, a region of the shocked solar wind between the 
termination shock (TS) and the heliopause (HP), has been ex- 
plored by Voyager 1 and 2 after the TS crossing. Voyager 1 
crossed the TS at a distance of 94 AU from the Sun in December 
2004 ( [Stone et al. |2005|l and Voyage r 2 crossed it at 84 AU in 
August 2007 (Richardson et al. 2008). The Voyager spacecraft 
revealed an important property of the outer heliosphere struc- 
ture - the asymmetry of the TS - and have been providing the 
first measurements of the heliosheath. Strong fluctuations of 
the magnetic field, the solar wind density, the velocity, and the 
temperature were observed by Voyager 2 in the heliosheath dur- 
ing about one year after the TS crossing (Richardson and Wang 
|2011| |Burl aga et al . |2010) . Fluctuation amplitudes have been 
decreasing as Voyager 2 moved deeper into the sheath. Richard- 
|son et aL| ( |2009| l also reported quasi-periodic variations of the 
solar wind velocity components and flow angles with a period of 
110 days observed by Voyager 2 in the heliosheath. The plasma 
experiment is not working on Voyager 1, but data from the mag- 



netic field experiment show large-scale magnetic field fluctua- 
tions in the heliosheath (Burlaga et al. 2006; Burlag almd Ness| 
f| |2010[ l. Observations from the Voyager spacecraft show that 
the heliosheath is a region with a highly variable and complex 
plasma flow. 

Different temporal and spatial effects may cause these dy- 
namic flows in the heliosheath. It is important therefore to in- 
vestigate the effects of temporal solar wind structures on the 
heliosheath to separate spatial and temporal plasma variations. 
Several solar wind phenomena contribute to temporal variations: 
the 11 -year solar cycle variations, large-scale structures such as 
interplanetary coronal mass ejections (ICMEs) (Richard son et| 
|aL]|2006a ; Richards on et al. l|2005|l and merged interaction re- 



gions (MIRs) (JBurlaga et al. |1993| l that usually occur at solar 



maximum; corotating interaction regions (CIRs) formed at the 
declining phase of solar activity (Burlaga et al. ||1997 20Q3|>; 
and CME- and CIR-driven shocks dBurlaga 1 1994) [Wang et al' | 



2001 |Wang and Richardson ||2002| |Richardson et al. ||20 06b). 



These structures generate significant changes in the solar wind 
parameters. Measurements show that during an 1 1 -year solar cy- 
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cle the solar wind ram pressure changes by factor of 2 from solar 
minimum to maximum. CIRs are characterized by an enhanced 
magnetic field, plasma density, and pressure, and are bounded 
by a pair of shocks. Observations of CIRs by Voyager 2 and Pi- 
oneer 10 in the inner heliosphere within 10 AU showed that the 



ram pressure in CIRs may increase by factor of 15-30 (Gazis 
2000 ). At larger heliospheric distances, CIRs expand and merge 



to form CMIRs ( Burlaga]T983| |Burlaga et aL 2003 1. Voyager 2 



data during the solar minimum in 1994-1995 near 45 AU show 
the sequence of recurring sharp, shock-like increases in the solar 
wind speed that resemble very much forward shocks (Lazarus et 



al. 1 199 9 ). These structures are associated with CMIRs. Changes 



of the solar wind parameters were smaller, but the ram pressure 
changed by factor of 2-4 at the shocks. Voyager 2 observations 
of solar maximum plasma between 65-70 AU indicated peri- 
odic fluctuations with a correlated solar wind speed, density, and 
magnetic field, which increases the solar wind ram pressure by a 
factor of 10 within timescales of 6-12 months (Richardson et al. 
[] ( |2QQ3[ >). Supposedly, these structures are candidates for MIRs, 
which are formed from merging transient solar wind events. All 
these solar wind structures propagate to the TS and beyond and 
affect the solar wind flow in the heliosheath. 

Effects of the 11 -year solar cycle variations on the interac- 
tion of the solar wind with the local interstellar medium (LISM) 



were studied by many authors (Karmesin et al. |1995 Bara 



and Malama 
Muller 12003 



2004a|bt 



nov and Zaitsev |19 98 Sche rer and Fahr ||2003a|b[ |Izmodenov 



Izmodenov et al. |2005||2b081|Zank and 



Pogorelov et al. |2009| l. Models applied different 



boundary conditions that simulated changes of the solar wind 
dynamic pressure during the solar cycle. Some global models 
that considered the solar cycle effects ignored the interstellar H 
atoms or used simplified fluid or multi-fluid approximations to 
describe the neutrals. A kinetic description for interstellar H 
atoms (which is necessary because of their large mean free path) 
was employed in Izmo denov et al. | ([2005 2008). The models 
showed that the TS and HP oscillate in response to varying solar 
wind ram pressure in the solar cycle. Using realistic boundary 
conditions in a time-dependent 2D kinetic-hydrodynamic model, 
Izmodeno v et al. | p008) obtained that the TS reflects variations 
of the solar wind dynamic pressure in 1-1.2 year and the TS lo- 
cation fluctuates by ± 7.5 AU. The model with 1 1-year periodic 
boundary conditions by Izmodenov et al. ( 2005 i showed that 
the HP varies less - by + 2 AU near its mean value. Therefore, 
the boundaries of the heliosheath are constantly in motion, and 
the heliosheath is expected to be a highly dynamic region. 

Interaction of the TS with various interplanetary distur- 
bances from upstream was studie d by|Barnes |(|1993[l,|Naidu and] 



Barnes ( 1994a ]b|, |Steinolfson |(fl994|l,|Story and Zank l ( |1995[ 
1997|l,|Baranov et al.|(|1996a|b|),|Ratkiewicz et al. | ( |1996| ), |Zank| 
and Muller | ( |2003| l, |Baranov and Pushkar | ( |2004| l and others in 
ID and 2D hydrodynamic and MHD models. The models stud- 
ied the propagation of solar wind shock waves, contact discon- 
tinuities, forward-reverse shock pairs, and ram pressure pulses 
through the TS and predicted the flow structure downstream of 
the TS. Recent works by |Washimi et al. | ( [2007] [20TT) explored 
the effects of realistic pulses of the solar wind ram pressure on 
the heliosheath flow. They performed a 3D MHD simulation us- 
ing Voyager 2 plasma data for the boundary conditions upstream 
of the TS. They showed that when the ram pressure pulse col- 
lides with the TS, (1) the TS moves away from the Sun; (2) a 
large-amplitude magnetosonic wave is generated downstream of 
the TS; (3) the magnetosonic wave is reflected inside the he- 
liosheath; and (4) the collision of the reflected wave with the TS 
causes the motion of the TS toward the Sun. 



Effects of CIRs in the heliosphere were modeled by |Pizzo| 
and Gosling"] ( |T994] l, |Pizzo | ( |1994a|b[) in a 3D MHD coro tating 
model of the solar wind flow. IBorovikov et al/H ( |2012| l mod- 
eled the evolution of CIRs in the outer heliosphere and in the 
heliosheath in a 3D time-dependent model of the solar wind in- 
teraction with the LISM. Their results show that CIRs create a 
complex non-stationary plasma flow in the heliosheath and pro 
duce entropy and fast magnetosonic waves. Burlaga et al. 



( |2Q03) > investigated the evolution of the magnetic field fluctua- 
tions induced by CIRs in a ID MHD model using realistic solar 
wind parameters in 1995 during the solar minimum. The model 
predicted broad regions with enhanced magnetic field caused by 
CIRs in the outer heliosphere. 

MIRs and global MIRs (GMIRs) are non-periodic large- 
scale disturbances in the heliosphere compared to CIRs. MIRs 
are regions with enhancements in the density and magnetic field 
strength and an increase of the bulk speed (Bur laga et al. |1993~ 



Burlaga a nd Ness"fl 994). These structures usually evolve from 
transient events that occur near the Sun, such as CMEs and iso- 



lated fast streams ( 


Richardson et al. |2006b Whang etal. 2001 


Richardson et al. 


2002). Interplanetary CMEs and MIRs were 



frequently observed by the Wind, Ulysses, and Voyager space- 
craft. Observations from Ulysses showed that these structures 
may drive a pair of forward and reverse shocks ( Whang et al. | 



[20UT] |Gosling et al. |1994| |Manchest e~and Zurbuchen |2006| ). 
Before crossing the TS, Voyager 2 observed a MIR at a distance 
of 79 AU fro m the Sun ([Richardson et al.||2006a| ). Using a ID 
MHD model, [Richardson et al.| ( |2006a| ) showed that the MIR was 
produced by a CME merged with high-speed plasma streams. 
The model also showed that the MIR is bounded by forward and 
reverse shocks which agrees with Voyager 2 observations. 

In general, the structure of a MIR is very complex because 
of interaction with ambient plasma and merging with other dis- 
turbances in the solar wind. In this work, we aim to explore 
the dynamical effects of a pair of shocks relevant to MIRs that 
propagate to the outer heliosphere and heliosheath. We describe 
in detail the interaction of a pair of shocks with the TS, prop- 
agation in the heliosheath, and the interaction with the HP and 
explain the formation of MHD discontinuities and waves in the 
heliosheath that create the highly variable flow. We use a 3D 
MHD global model of the interaction between the solar wind 
and the LISM ( |Opher et al. |2009] ). 

Previously, interaction of a pair of forward-reverse shocks 



with the TS was studied by|Story and Zank (1995 |1997)inthe 



Baranov 



frame of ID planar gas-dynamic and MHD models. 
[eTaLl ( |1996a|b| ), |Baranov and Pushkar | p004[ ) considered 2D 
MHD interaction of the TS with forward and reverse interplan- 
etary shocks. They showed that the collision of shocks with the 
TS generates a large number of wave modes and shocks down- 
stream of the TS. We compare our results with the previous mod- 
els and extend the study to the interaction of shocks with the HP. 

The outline of the paper is the following. In section 2 we 
briefly describe the model. In section 3 we present the results: 
(1) the evolution of a pair of forward-reverse shocks in the su- 
personic solar wind, (2) the interaction with the TS, and (3) the 
propagation through the heliosheath and interaction with the HP. 
Conclusions are given in section 4. 



2. Description of the model 

To model the propagation of the solar wind disturbances to the 
heliosheath we used a global 3D MHD model of the solar wind 
interaction with the LISM ( {Opher et a"l~| |2"009). The model is 
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based on the BATS-R-US code ( |Toth et al. |20l2) . The model in- 
cludes the influence of interstellar hydrogen atoms that penetrate 
into the heliosphere. Plasma protons interact with H atoms in a 
charge-exchange process. Because of the large mean free path 
of H atoms compared to the size of the heliospheric interface, a 



B: 0.050 0.105 0.161 0.216 0.271 0.327 0.382 0.438 0.493 0.548 



kinetic approach is required for the neutral component ( Izmode- 
|novetal. |2000||rzmodenov~p 001 ). Several global models use a 
multi-fluid approximation to describe H atoms ( |Zank andM uller 
2003] Opher et al. |2009 1 due to the computational difficulty of 



solving a kinetic equation. The multi-fluid model reproduces the 
global structure of the outer heliosphere well, but it has some 
limitations (Alexashov and Izmod enov |2005[ ). The comparison 
of the multi-fluid and kinetic models performed by | Alexashov] 
and Izmodenov ( 2005) > showed that solutions for the plasma dif- 



fer by 5%. The difference in the parameters of the interstellar 
hydrogen atoms between the kinetic and multi-fluid models is 
significant. In the present work, we are interested in the plasma 
flow and used the multi-fluid approximation for the description 
of H atoms. A study of non-stationary flows in the heliosheath 
based on a 3D MHD model with a kinetic treatment of H atoms 
(as was done by Alouani-Bibi et al. ( 201 ![ )) will be performed 
in the future. 

Separate systems of Euler equations with source terms for 
charge exchange were solved for each of the four H atom popu- 
lations, and ideal MHD equations with source terms were solved 



for the plasma component (for more details, see Opher et al. 



(2009 )). The solar and interstellar magnetic fields are included in 
the model. The inner boundary of a domain is a sphere at 30 AU 
and the outer boundary is at x = + 1000 AU,y - ±1000AU,z = 
+ 1000AU. Parameters of the solar wind at the inner bound- 



ary we re chosen to match the values obtained by Izmodenov 
(2009) at 30 AU: V sw = 417km/ s , « - c ™ - 



= 8.74 x \0-=cm- 



T sw = 1.087 x 10 K and the Parker spiral magnetic field B sw = 
1 Al X 10" 3 «r at the equator. Our simulations are different from 



Opher et al. ( 2011| l in that we assumed that the magnetic axis 
is aligned with the solar rotation axis. For a steady-state solu- 
tion of the interaction between the solar wind and the LISM, 
the solar wind flow at the inner boundary was assumed to be 
spherically symmetric. For the interstellar plasma we assumed: 
Vlism = 26Akm/s, n U sM = 0.06cm -3 , and Tlism = 6519K. 
The number density of H atoms in the interstellar medium is 
n n U sM - 0.18cm -3 , the velocity and temperature are the same 
as for the interstellar plasma. The coordinate system is such 
that the Z-axis is parallel to the solar rotation axis, the X-axis 
is 5° above the direction of interstellar flow and Y completes the 
right-handed coordinate system (a schematic figure can be found 
in |Alouani-Bibi et al. | ( |20TT] l). 

We assumed that the interstellar magnetic field has a magni- 
tude Busm - 4.4/iG and the orientation is such that the angle 
between the interstellar flow velocity Vlism and Busm is 20°, 
and the angle between the (Bus m, Vlis m) plane and the solar 
equator plane is 60°. This direction for the interstellar magnetic 
field is close to the hydrogen deflection plane ( |Lallement et al.| 
2005 ) and reproduces the TS asymmetry in the Voyager 1 and 
Voyager 2 directions ( |Opher et al. |2009 1. 



In this paper, we consider the propagation of the solar wind 
disturbances along the Voyager 2 trajectory. We designed a nu- 
merical grid with a highly refined cone extending from 30 AU 
beyond the HP along the Voyager 2 trajectory. The entire grid 
in the domain has 10 7 cells with nine levels of refinement, rang- 
ing from scales of 0.5 AU at the inner boundary and in the cone 
to 32 AU at the outer boundary. To solve the equations numeri- 
cally, we used a second-order HLLE scheme with a monotonized 
central flux limiter function. 




Fig. 1. Plane cut through the solar rotation axis and Voyager 2 direc- 
tion (Oz-V2 plane) from a 3D MHD simulation showing the magnitude 
of the magnetic field (nT). The Voyager 1 trajectory is ~ 30° above 
the solar equator (white line) and that of Voyager 2 is ~ 30° below the 
equator (blue line). Black boundaries show the regions with different 
levels of grid refinement, the resolution in the cone along the Voyager 2 
trajectory is 0.5 AU. 



Fig. [T]shows the magnitude of the magnetic field in the plane 
through the solar rotation axis and the direction of the Voyager 2 
trajectory (hereafter Oz-V2 plane) for the steady-state solution. 
The inclined Busm produces an asymmetry of the heliosphere 
that pushes the TS and HP closer to the Sun in the south. This 
stationary solution was used to initiate the disturbances at the 
inner boundary. 

To generate a pair of forward and reverse shocks we in- 
creased the solar wind speed from 4\lkmjs by factor of a 1.5 
at the inner boundary along the Voyager 2 direction. We as- 
sumed that the angle between the discontinuity surface and the 
solar wind velocity vector in the direction of Voyager 2 is 90 
degrees (see Fig. 2 below). The angular size of a high-speed 
stream is 30° in 9 and 30° in (p (where and <p are the latitudinal 
and azimuthal angles in the HGI coordinate system). The shocks 
and the interaction region between them are located within the 
highly resolved grid cone during the entire time of the simula- 
tion. We split the study of the propagation of the pair of forward- 
reverse shocks into the following sections: 1) propagation in the 
supersonic solar wind out to the TS; 2) interaction with the TS; 
3) propagation in the heliosheath; and 4) interaction with the 
HP and the posterior evolution of the reflected waves in the he- 
liosheath. 



3. Results 

3.1. Evolution of a pair of shocks in the supersonic solar wind 

After the initiation of a high-speed stream in the Voyager 2 di- 
rection, an arbitrary discontinuity separating the fast and am- 
bient slow solar wind forms at the front. It decays and two 
shocks - forward and reverse - are formed. Figure [2] presents 
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log,„n,cm 3 ; -2.7-2.6-2.4-2.3-2.2-2.1-1.9 V, km/s: 50 124 198 271 345 419 493 567 




log lc p: -13.0-12.7-12.5-12.2-11.9-11.6 B, nT: 0.05 0.13 0.22 0.30 0.38 0.47 0.55 




Fig. 2. Cut from a 3D MHD simulation showing the log of the solar 
wind number density (em -3 ), the velocity (km/s), the log of thermal 
pressure (dyn/cm 2 ), and the magnitude of the magnetic field (nT) in 
the Oz-V2 plane. The plots show the formation of a forward-reverse 
shock pair that propagates along the Voyager 2 direction (dashed white 
line). The forward shock (FS), reverse shock (RS), and termination 
shock (TS) are denoted. Black boundaries show the regions of the grid 
refinement. 

the solar wind number density, the velocity, the thermal pres- 
sure and the magnetic field in the Oz-V2 plane at t — 0.3 yr 
(f = corresponds to the steady-state solution at the time that 
the disturbance was initiated). It can be seen that a large-scale 
interaction region with enhanced density, thermal pressure, and 
magnetic field is formed in the supersonic solar wind; the region 
is bounded by the forward and reverse shocks. Indeed, a solu- 
tion of the Riemann problem in MHD shows that a tangential 
discontinuity (TD) must form between the two shocks. In fig- 
ure 3 (blue curve) an increase of plasma density and magnetic 
field and decrease of temperature can be identified between the 
shocks. These parameter changes correspond to the TD. The 
numerical scheme used in the simulation diffuses tangential dis- 
continuities. For this reason we focus on the evolution of the two 
shocks in our MHD analysis of the plasma flow. 

At t — 0.3 yr the interaction region passes the distance 45 
AU from the Sun and the compression ratios are 5fs - 2.4 
and 6rs = 2.1 for the forward and reverse shocks, respectively. 
The total pressure p, , = p,/, + pg and dynamic pressure pV 2 
increase by factors of 5 and 4, respectively, in the interaction 
region compared to the undisturbed solar wind upstream of the 
forward shock. 

The characteristics of the shocks and the variations of the 
solar wind parameters in the interaction region change while 
they propagate away from the Sun. The shocks weaken due to 
a change of the background solar wind parameters with helio- 
spheric distance and spherical expansion of the interaction re- 
gion. Our model shows that the forward shock compression ra- 
tio 6fs decreases by 13 % from 40 AU out to the TS. For the 
reverse shock, 6rs decreases by ~ 20%. The forward shock 
speed decreases by ~ 60% from 288 to 120 km/s near the TS 
(with respect to the solar wind upstream of the forward shock). 



3 Or °- 39yr 1 - 6 r 

0.42 yr 




Fig. 3. Profiles of normalized density p/po, velocity V/V , temperature 
T/Tq, and magnetic field B/Bq along the Voyager 2 trajectory for dif- 
ferent times showing the change in amplitudes of the fluctuations in the 
interaction region. Here po, Vo, Tq Bq denote the steady-state solution. 
Notations: FS - forward shock, TD is the tangential discontinuity, RS - 
reverse shock. 



The reverse shock accelerates in the solar wind frame, the shock 
speed increases from 120 km/s to 245 km/s. 

The interaction region expands radially with time since the 
forward and reverse shocks move in opposite directions in the 
solar wind frame. The width increases about three times from 
4.5 AU at a distance ~ 40 AU from the Sun to 13.5 AU near the 
TS. 

To explore the 3D effects of the flow, we analyzed the vari- 
ation changes of the solar wind parameters in the interaction re- 
gion while it propagates in the supersonic solar wind. To reveal 
the latitudinal dependence we investigated the different direc- 
tions in the V2-Oz plane: 6\ = 30° corresponding to the Voy- 
ager 2 direction; 62 = 20°; #3 = 40° (here the angle is measured 
from the equatorial plane toward the south). These directions are 
within the high resolved cone of our computational grid. Figure 
[3] shows the evolution of the interaction region along the Voy- 
ager 2 direction. Plasma parameters vary inside the interaction 
region because of the TD and possibly other waves. A region of 
highest density and magnetic field is created through the pile-up 
of the solar wind plasma behind the forward shock. The inter- 
action region between the shocks exhibits strong variations of 
the solar wind parameters. At a distance of 45 AU from the 
Sun, the magnetic filed strength increases by a factor of 2.4 in- 
side the interaction region, the solar wind density by 2.4, and the 
temperature by 2.5. As the disturbance evolves, the fluctuation 
amplitude of normalized magnetic field B/Bo increases by 6%, 
the amplitude of the density fluctuation p/po decreases by 12%, 
the temperature fluctuation 77 To increases non-monotonically 
by 10%, and the change in the velocity V/Vo is negligibly small. 
Here po, Vo, To and Bq refer to the steady-state values. Along 
the directions 62 and #3 the model shows the same fluctuations 
behavior for p/po, V/Vo, and p to tlptoi - However, in the direc- 
tion 63 the amplitude of B/B^ fluctuation increases by 8%, while 
along the 82 it remains nearly constant. The difference in evolu- 
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tion of the magnetic field fluctuations is caused by the variation 
of the background Parker magnetic field Bo with latitude. Lati- 
tudinal variations in evolution of the interaction region show that 
the 3D effects of the flow induced by realistic large-scale solar 
wind disturbances could be important. 

Our model shows increasing magnetic field fluctuations in 
the interaction region with distance, which is different from what 
orlaga et al. (2003 1 found in their ID model. They found 
that the amplitude of B/Bq in CIRs decreases by ~ 4%. This 
discrepancy is due to a difference in the boundary conditions 
between the models and 3D effects taken into account in our 
model. 

Modeling the evolution of an interaction region bounded by 
a pair of shocks from 30 AU out to the TS showed that (1) the 
shocks weaken; (2) the interaction region expands radially and 
decelerates; (3) the variations of the solar wind parameters in 
the interaction region become weaker, and the variations of the 
magnetic field increase with the heliospheric distance. Our study 
of the 3D effects of the solar wind flow due to the propagating 
pair of shocks showed that magnetic field fluctuations become 
stronger with increasing latitude from the Voyager 2 direction. 

3.2. Interaction of a pair of shocks with the TS 

When a pair of forward-reverse shocks propagates to the outer 
heliosphere, the forward shock eventually encounters the TS. 
Bara nov et al . (1996a) studied MHD interaction of a for- 
ward interplanetary shock with the stationary TS in a two- 
dimensional model. They showed that the solution of the prob- 
lem is determined by five dimensionless parameters: Mps = 
Vrs/ao, Mip S = V/psM), B - %npjB 2 ,Q, and ifr TS , where Vjs 
is the solar wind speed immediately upstream of the TS, Oq is 
the sound speed upstream of the TS, Vjps is the interplanetary 
shock speed (in our case this is the speed of the forward shock), 
B is the ratio of thermal and magnetic pressure, 8 is the angle be- 
tween the normals to the TS and the interplanetary shock fronts, 
and i//rs is the inclination angle of the interplanetary magnetic 
field vector relative to the TS front. 

In our simulation, 9 ~ 0, if/js ~ 0, Mrs = 3.8, Mjps = 3.5, 
B = 4.6. For these values, |Baranov et al.| ( | 1996a} reported a 
configuration of the interaction of the forward shock and the TS 
that was expressed by the following scheme: 



log 10 n, cm" 3 : - 2.8 -2.7 -2.5 -2.4 -2.3 -2.2 -2.1 V, k m/s: 50 125 200 275 350 425 500 575 
•20 - 



TS 



FS 



FS' <- R C SS TS' 



where the arrows show the directions of motion of the discon- 
tinuities after the shocks interact in the solar wind frame. Here 
FS' is the new modified forward shock, propagating into the 
heliosheath, R is a slow rarefaction wave, C is a contact discon- 
tinuity, SS is a reverse slow shock, and TS' is the modified TS. 

Our numerical solution shows structures in the flow similar 
to the analytical solution of Baranov et al. (1996a ). Figure |4] 
shows the number density, the velocity, the temperature, and the 
magnetic field in the V2-Oz plane at the moment t = Q.lyr after 
the forward shock crossed the TS and before the reverse shock 
encounters the TS. Our results show the formation of a new for- 
ward shock (FS') that propagates into the heliosheath, the modi- 
fied TS (TS'), and a structure between them characterized by an 
increasing plasma density and magnetic field intensity accompa- 
nied with a decreasing temperature (and thermal pressure) and 
no change in velocity. At the same time, the total pressure is 
constant across this structure, indicating that this is a tangential 
discontinuity (denoted as TD in Fig. 4). 

Sams onov et al. | (p006 ) considered an interaction of an in- 
terplanetary shock with the Earth's bow shock. Their solution 




-6° - 40 -80 -60 -40 

, J, K: 13.4 13.8 14.1 14.4 14.7 15.1 15.4 B, nT: 0.08 0.12 0.17 0.24 0.31 0.38 0.44 0.51 
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Fig. 4. Meridional cuts through the Oz-V2 plane from a 3D MHD 
simulation showing the log of the solar wind number density (cm~ 3 ), 
the velocity (km/s), the log of temperature (K), and the magnetic field 
magnitude (nT) at t = 0.7 yr when the forward shock crossed the helio- 
spheric termination shock. FS' denotes the forward shock propagating 
into the heliosheath, TD is a tangential discontinuity, TS is the termina- 
tion shock, and RS is a reverse shock upstream of the TS. 



shows the same qualitative behavior of the plasma parameters 
between the two shocks after their interaction. The authors sug- 
gested that this structure is a combination of a slow rarefaction 
wave, a contact discontinuity, and a slow shock. The increas- 
ing magnetic field with the decreasing density is generated by 
the reverse slow shock. This decreasing density is compensated 
by a stronger increase of density across the contact discontinu- 

(2006) 



ity. Velocity variations are negligible. Samsonov et al. 



pointed out that the increased grid resolution in their model is 
not sufficient to reproduce the separated discontinuities because 
of similar velocities of the discontinuities and small changes of 
MHD parameters. 

Even with the second-order numerical scheme and the spa- 
tial resolution of 0.5 AU, our simulation is unable to resolve such 
small structures. If the slow shock exists in our case according to 
the solution of Baranov et al. (1996a]), it may not be resolved by 
our numerical method and the grid used. High-resolution runs 
were performed by Op her et al. | ( |201 1\ , but they are computa- 
tionally extremely costly. 

Since the TS is a reverse shock in the solar wind frame of 
reference, its interaction with the incident forward shock results 
in a weakening of both shocks. In our case the strength of the 
forward shock decreased by ~ 30% after the TS crossing. The 
strength of the TS decreased by 7%. Due to the increasing ram 
pressure behind the forward shock, the TS is displaced by 3 AU 
away from the Sun. 

After the passage of forward shock the modified TS' inter- 
acts with the reverse shock - the rear side of the interaction re- 
gion. MHD interaction of the TS with the reverse interplanetary 
shock was considered by Baranov and Pushkar (2004 ) in a two- 
dimensional model. Their solution gives a configuration consist- 
ing of a new TS and several waves propagating downstream of 
the TS - a fast rarefaction wave, a slow shock, a contact disconti- 
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Fig. 5. Meridional cuts through the Oz-V2 plane from a 3D MHD 
simulation showing the log of the solar wind number density (cm -3 ), 
the velocity (km/s), the log of the temperature (K), and the magnetic 
field magnitude (11T) at t = 1 yr after the reverse shock crossed the TS. 
FS' denotes the forward shock propagating into the heliosheath; TD\ 
and TD 2 are tangential discontinuities; TS is the termination shock; R 
denotes a rarefaction wave formed after the interaction of the TS and 
the reverse shock. 



nuity, and a slow rarefaction wave. Figure[5]shows the results of 
our simulation; the solar wind parameters are plotted as in Fig- 
ure |4]but at t — lyr when the reverse shock crossed the TS'. It 
can be seen that a tangential discontinuity (denoted TD2) and a 
rarefaction wave R propagate into the heliosheath after the TS' 
interaction with the reverse shock (Fig. 5). At the tangential dis- 
continuity, the plasma density and the magnetic field decrease 
and the temperature increases, while the total pressure remain 
constant. In the rarefaction wave R, the parameters change con- 
tinuously. The velocity change in the rarefaction wave is negli- 
gibly small. Since the configuration of a shock-shock interaction 
in MHD could be very complex, some additional waves and dis- 
continuities may arise in the heliosheath that are not resolved in 
the simulation. Our simulation shows that after merging with the 
reverse shock, the strength of the TS increases by 10% and the 
TS is displaced by 2.6 AU away from the Sun due to increasing 
solar wind dynamic pressure at the reverse shock. 

Therefore, when an interaction region bounded by a pair of 
shocks passes the TS, a new large-scale region of disturbed solar 
wind plasma forms in the heliosheath and propagates toward the 
HR This new disturbance has a complex structure, it is bounded 
by the forward shock FS' at the front and the tangential discon- 
tinuity TD2 at the rear (see Fig. [5]). The TS strength is little 
changed by the passage of this pair of shocks. The forward 
shock weakens after interacting with the TS. In the heliosheath, 
the forward shock strength varies with latitude: the values of 
6fs obtained within ±10° from the Voyager 2 direction showed 
that Sfs remains nearly constant toward the higher latitudes and 
decreases by 5% to the equator. At the tangential discontinuity 
TD2, the solar wind number density decreased by 30%. The re- 
gion between FS' and TD2 is compressed heated plasma with 
an enhanced density and magnetic field, moving into the he- 



Fig. 6. Meridional cuts in the Oz-V2 plane from a 3D MHD simulation 
showing the distribution of the solar wind density (cm -3 ) for different 
times when the interaction region propagates in the heliosheath. The 
plot (a) corresponds to the moment of time t = lyr; (b) t = 1.1 yr; (c) 
t= 1.2yr;and(d)f = 1.3 yr 



liosheath faster than the ambient plasma (plot of velocity in Fig. 
[5}. The solar wind is highly variable in this region due to for- 
mation of waves and discontinuities, as shown above. In the 
following section, the propagation of this large-scale disturbed 
region through the heliosheath is discussed. 



3.3. Propagation in the heliosheath 

The heliosheath is a region with spatial variations of the solar 
wind parameters that may affect the shock properties and ampli- 
tudes of the disturbances propagating from the supersonic solar 
wind. Along the Voyager 2 direction from the TS to the HP, the 
plasma decelerates and the magnetic field magnitude increases. 
A region with enhanced magnetic filed exists near the HP (see 
Fig. [TJ. In our simulation the forward shock transmitted from 
the supersonic solar wind into the heliosheath is a weak perpen- 
dicular shock. For a weak perpendicular shock a dependence of 
the shock speed on the background parameters is given by (Gur- 
|nett and Bh attacharj ee |2005) 



V 2 = 2(c 2 + c 2 )f(8- 1)(<W~1), 

where V„ = (V sw - D) • n is the normal component of upstream 
plasma velocity in the shock frame, D is the shock speed, c s 
is the sound speed upstream the shock, c A is the Alfven speed 
upstream the shock, 6 is the ratio of upstream and downstream 
densities at the shock, and § mclx - (y + l)/(y- 1), where y — 5/3. 
In the heliosheath along Voyager 2, c s decreases and c A in- 
creases due to increase of magnetic field. The term (c 2 + c\) 
remains nearly constant across the heliosheath, which gives a 
constant V 2 . But the solar wind velocity decreases across the 
heliosheath due to deceleration toward the HP. Therefore, a de- 
celeration of the forward shock in the heliosheath should occur. 
Our simulation shows that the forward shock speed decreases by 
10%. While the large-scale disturbance moves through the he- 
liosheath, the compression of the shock also decreases by 6%. 
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Fig. 7. Profiles of the normalized density p/po, the total pressure 
Piot/Ptoio, the temperature T/T , and the magnetic field B/B along the 
Voyager 2 trajectory for different moments of time showing the evolu- 
tion of the fluctuations when the interaction region propagates in the 
heliosheath. Notations: FS - forward shock, TS - termination shock, 
TD U TD 2 - tangential discontinuities. 



The 3D effects of the flow become more pronounced when the 
forward shock approaches the HP. Figure [6] shows plasma den- 
sity and velocity streamlines in the V2-Oz plane in different mo- 
ments of time while the disturbance propagates through the he- 
liosheath. Near the HP 6fs = 1.3 in Voyager 2 direction. Values 
of Sfs within ±10° from the Voyager 2 direction show a latitu- 
dinal dependence: 6fs increases by 5% with increasing latitude 
and decreases by 10% from the Voyager 2 toward the equator. 
The forward shock becomes stronger with increasing latitude, 
because the solar wind turns around the heliopause to the he- 
liospheric tails and compressed plasma behind the shock moves 
to the higher latitudes (Fig. [§}. Our analysis of 3D effects was 
performed within the narrow highly resolved grid cone, but in 
general stronger 3D variations may exist in the heliosheath. 

Our simulation shows an essential radial and latitudinal ex- 
pansion of the interaction region in the heliosheath. After the TS 
crossing, the radial width of the interaction region is 22 AU; it 
increases by 60% closer to the HP. From Fig. [6] it can be seen 
that the disturbed region significantly spreads in latitudinal ex- 
tent. The interaction region with latitudinal extension ~ 30° in 
the supersonic solar wind produces a large-scale disturbed region 
in the heliosheath with an extension of about 50°. 

Figure [7] shows the normalized solar wind number density, 
the total pressure, the temperature, and the magnetic field along 
the Voyager 2 direction for different moments of time while a 
disturbance travels throughout the heliosheath. The disturbance 
is characterized by two peaks of the magnetic field: B increases 
at the forward shock FS and at the tangential discontinuity TD\ . 
At the forward shock the increase of B/Bq is about 30%, at 
the tangential discontinuity the magnetic field increases stronger, 
~ 60%. Density profiles show that the disturbance also creates 
two regions with maximum plasma density. The amplitudes of 
p/po in both regions are comparable, ~ 30%. As the disturbed 
region moves toward the HP, the amplitude of B/Bq fluctuation 




Fig. 8. Meridional cuts through the the Oz-V2 plane from a 3D MHD 
simulation showing the contours of SB/B for different times. Nota- 
tions: FS is the forward shock, HP is the heliopause, FW is reflected 
fast magnetosonic wave propagating toward the TS, and TS is the he- 
liospheric termination shock. The temporal evolution is shown in an 
animation available in the online edition. 



decreases by 12 %, p/po decreases by 1 1 %, and V/Vq increases 
by 20 % along the Voyager 2 direction. 

One year after crossing the TS, the forward shock encounters 
the HP (left top plot in Figure 8). The interaction of the forward 
shock with the HP causes a motion of the HP away from the 
Sun by few AUs. As the result of the shock interaction with 
the HP, a new shock is created that propagates into the inter- 
stellar medium, a magnetosonic wave is reflected from the HP 
inside the heliosheath and propagates back toward the TS. Fig. [8] 
presents the propagation of the reflected wave in the heliosheath 
toward the TS. It shows contours of SB/Bq = (B - Bo) / Bo, where 
Bo refers to the steady-state solution. FS denotes the forward 
shock that propagates in the interstellar medium. FW denotes 
the magnetosonic wave reflected from the HP. The locations of 
TS and HP are shown in the left top plot in Fig. 8. The re- 
flected magnetosonic wave FW is a fast magnetosonic wave. 
From our simulation, the average fast magnetosonic speed in the 
heliosheath is about 250 km/s. Using the results presented in 
Figure 8, one can calculate that the speed of reflected magne- 
tosonic wave is about 230 km/s. This speed is comparable to the 
fast magnetosonic speed in the heliosheath. It can also be seen 
that vortices may form in the heliosheath at the sides of the high- 
speed stream. The formation of vortex flow and possible role of 
K-H instability will be discussed in subsequent paper. 

The reflected magnetosonic wave FW approaches the TS. 
The interaction of the TS with FW occurs about 2.5 years af- 
ter the forward shock entered the heliosheath. A bounce of the 
magnetosonic wave causes the displacement of the TS by 2-3 AU 
toward the Sun. Because of the complex flow in the heliosheath 
it is hard to identify the secondary reflected waves from the TS. 
However, the formation of a cascade of small disturbances in so- 
lar wind plasma reflecting from the TS and propagating back 
to the heliosheath is seen from the simulation (see animation 
linked to Figure |8]in the on-line material). After the interaction 
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with the magnetosonic wave the TS moves into the heliosheath. 
This movement may also contribute to the formation of com- 
pressional waves in the solar wind plasma downstream of the 
TS, since the TS acts like a piston that pushes the solar wind 
material inside the heliosheath. These secondary magnetosonic 
waves reflected from the TS are caught up by the bulk solar wind 
flow around the HP and are carried away to the tails of the helio- 
sphere. 

4. Conclusions 

Using a global 3D MHD model of the solar wind interaction with 
the LISM, we studied a propagation of a pair of forward and re- 
verse shocks from the supersonic solar wind into the heliosheath. 
A pair of shocks is usually driven by MIRs that form in the solar 
wind during solar maximum. The pair of shocks was initiated 
in the direction of the Voyager 2 trajectory. To capture newly 
created discontinuities formed in the solar wind through the pas- 
sage of a pair of forward-reverse shocks, we used a non-uniform 
spatial grid with highest resolution of 0.5 AU in the cone along 
the Voyager 2 trajectory. 

To study the propagation of the pair of shocks from the su- 
personic solar wind to the HP we analyzed (1) the evolution in 
the supersonic solar wind, (2) the interaction with the TS, (3) the 
propagation through the heliosheath, and (4) the interaction with 
the HP. The pair of shocks is formed by a steep increase of the 
solar wind speed at 30 AU from 417 km/s to 625 km/s. The sim- 
ulation in the supersonic solar wind showed that the interaction 
region between the shocks expands and decelerates; both shocks 
weaken while moving to the larger heliospheric distances. When 
the interaction region is at 45 AU from the Sun, the plasma den- 
sity p, the temperature T, and the magnetic field strength B vary 
inside the region by factors of 2.4, 2.5, and 2.3. While the in- 
teraction region propagates outward, the variation amplitudes of 
the plasma parameters change: in the Voyager 2 direction the 
density fluctuations p/po decrease by 12 %, T/Tq increases by 
10 %, B/Bq increases by 6 %. We found that magnetic field 
fluctuations behave differently depending on the latitude. 

Modeling the propagation of a shock pair in the heliosheath 
showed the following effects: the collision of the pair of forward 
and reverse shocks with the TS causes the motion of the TS away 
from the Sun of a few AU and several new discontinuities are 
created downstream of the TS: a fast forward shock, tangential 
discontinuities, and possibly rarefaction waves and slow shocks. 
The passage of the forward shock through the TS results in a 
weakening of both shocks. The structure of the interaction re- 
gion between the shocks changes after the TS crossing: in the 
heliosheath it is bounded by a weak forward shock at the front 
and a tangential discontinuity at the rear. Newly formed dis- 
continuities in the heliosheath cause variations in the solar wind 
parameters and disturb the heliosheath flow. The magnetic field 
strength in the interaction region increases by 30-60%, the den- 
sity by 30%, the temperature by 30%. 

While the forward shock propagates toward the HP, the 
shock strength decreases and the shock decelerates. Magnetic 
field and density fluctuations have smaller amplitudes as the in- 
teraction region propagates deeper in the heliosheath. The in- 
teraction of the forward shock with the HP causes the outward 
motion of the HP and reflection of fast magnetosonic wave from 
the HP. Reflected waves propagate inside the heliosheath toward 
the TS and encounter the TS after about 1 .5 year. The interaction 
of magnetosonic waves with the TS displaces the TS toward the 
Sun and generates a secondary reflection of the magnetosonic 
waves from the TS. The propagation of reflected waves between 



the TS and the HP contributes to the dynamic flow in the he- 
liosheath solar wind. 

In this study we considered an evolution of a shock pair that 
was initialized by an abrupt increase of the solar wind speed. 
In a subsequent study we will focus on modeling realistic solar 
events that produce a pair of shocks or one strong interplanetary 
shock in the heliosphere driven by MIR. For example, series of 
solar events during August-September 2005 detected at the Earth 
produced a MIR and an associated strong shock that were ob- 
served at Voyager 2 immediately before the TS (at 79 AU) and 
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then at Voyager 1 beyond the TS. Webber et al. 



analyzed temporal variations of > 70MeV cosmic ray intensi- 
ties observed on Voyager 1 during the shock's travel through the 
TS and HP. Based on the data, they estimated the shock arrival 
times to the TS and the HP and the shock propagation speed in 
the heliosheath. A comparison of a realistic solar event simu- 
lation with observations would significantly improve our model 
and enable us to use it to predict heliosheath flows caused by the 
solar events during the present increase of solar activity. 
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